function g_k = numerical_gradient(sc, st, F, mu, W)

g_k = zeros(5,1);
I = eye(7);

h = 1e-9;

sc_p = sc(1);
sc_f = sc(2);
sc_g = sc(3);

sc_a = sc_p./(1-sc_f.*sc_f-sc_g.*sc_g);
sc_e = sqrt(sc_f.*sc_f+sc_g.*sc_g);

e = I(:,1);
sca = sc + h*e * sc_a*(1-sc_e^2);
scb = sc - h*e * sc_a*(1-sc_e^2);

Qa = proximity_quotient(sca, st, F, mu, W);
Qb = proximity_quotient(scb, st, F, mu, W);

g_k(1,1) = ( Qa - Qb ) / (2*h);


for i = 2:5
    
    e = I(:,i);
    sca = sc + h*e;
    scb = sc - h*e;
    
    Qa = proximity_quotient(sca, st, F, mu, W);
    Qb = proximity_quotient(scb, st, F, mu, W);
    
    g_k(i,1) = ( Qa - Qb ) / (2*h);
    
end



